
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <omp.h>
//////////////////////////////
#include <cuda.h>
#include <cuda_runtime.h>
//////////////////////////////

#define REPEAT 256

__device__ __inline__ float f0(float x,float y){
  return  sinf(x*x+y*y);
}

__device__ __inline__ float f1(float x,float y){
  return  cosf(x*x-y*y+1.0f/(x*x));
}

__device__ __inline__ float f2(float x,float y){
  return  tanf(x-y);
}

__device__ __inline__ float f3(float x,float y){
  return  log2f(x+y);
}

__device__ __inline__ float f4(float x,float y){
  return  logf(x*x+1.0f);
}

__device__ __inline__ float f5(float x,float y){
  return  log10f(x*x+2.0f/(y*y));
}

__device__ __inline__ float f6(float x,float y){
  return  sqrtf(1.0f+y*y);
}

__device__ __inline__ float f7(float x,float y){
  return  cbrtf(2.0f+y*y);
}


__global__ void non_divergent_kernel(float * __restrict__ a,
				     const float * __restrict__ b){
  int id   = blockIdx.x * blockDim.x + threadIdx.x;

  a[id] = f0(a[id],b[id]);
}

__global__ void divergent_kernel(float * __restrict__ a,
				 const float * __restrict__ b){
  int id   = blockIdx.x * blockDim.x + threadIdx.x;

  switch (threadIdx.x % 8) {
  case 0:
    a[id] = f0(a[id],b[id]);
    break;
  case 1:
    a[id] = f1(a[id],b[id]);
    break;
  case 2:
    a[id] = f2(a[id],b[id]);
    break;
  case 3:
    a[id] = f3(a[id],b[id]);
    break;
  case 4:
    a[id] = f4(a[id],b[id]);
    break;
  case 5:
    a[id] = f5(a[id],b[id]);
    break;
  case 6:
    a[id] = f6(a[id],b[id]);
    break;
  case 7:
    a[id] = f7(a[id],b[id]);
    break;
  }
}

int main(void){
  long unsigned int N = 1024L*1024*256;
  float *a = (float *)malloc(N*sizeof(float));
  float *b = (float *)malloc(N*sizeof(float));
  struct timespec ts0,ts1,ts2;
  double time;
  
  fprintf(stderr,"\nFilling arrays\n");
  for (long unsigned int i=0;i<N;i++){
    a[i] = (float)(rand()-RAND_MAX/2)/(float)(RAND_MAX);
    b[i] = (float)(rand()-RAND_MAX/2)/(float)(RAND_MAX);
  }

  
  // Allocate arrays on GPU
  float *d_a,*d_b;
  cudaMalloc(&d_a,N*sizeof(float));
  cudaMalloc(&d_b,N*sizeof(float));
  
  // Transfer initial data to GPU
  cudaMemcpy(d_a,a,N*sizeof(float),cudaMemcpyHostToDevice);
  cudaMemcpy(d_b,b,N*sizeof(float),cudaMemcpyHostToDevice);


  clock_gettime(CLOCK_REALTIME,&ts0);
  
  fprintf(stderr,"\tnon divergent code");
  for (int count=0;count<REPEAT;count++){
    non_divergent_kernel<<<N/128,128>>>(d_a,d_b);
  }
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-ts0.tv_sec + (ts1.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr," : %.2f\n",time);
 
  
  fprintf(stderr,"\tdivergent code    ");
  for (int count=0;count<REPEAT;count++){
    divergent_kernel<<<N/128,128>>>(d_a,d_b);
  }
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&ts2);
  time = (ts2.tv_sec-ts1.tv_sec + (ts2.tv_nsec-ts1.tv_nsec)/1e9);
  fprintf(stderr," : %.2f\n",time);

  exit(0);
}
